TP Integrador - Aprendizaje Estadístico

Enunciado

Consigna

En el archivo sat.txt se tiene datos de valores de espectros de pixels en una imagen de satélite, utilizados para la predicción de la clase de suelo. Realice un análisis de los datos utilizando diferentes métodos y compárelos mediante Validación Cruzada. Luego aplíquelos a la muestra de test sat.tst y analice los resultados obtenidos.

Propósito

La base de datos contiene los valores multi-espectrales de píxeles en regiones de 3x3 en una imágen satelital, y la clasificación asociada con el píxel central de cada región. El objetivo es predecir esta clasificación dados los valores multi-espectrales.

Detalle

Esta base de datos se generó a partir del escáner multi-espectral LANDSAT. Un marco de imágenes Landsat MSS consta de cuatro imágenes digitales de la misma escena en diferentes bandas espectrales. Dos de estos están en la región visible (correspondiente aproximadamente a las regiones verdes y rojas del espectro visible) y dos están en el infrarrojo. Cada píxel es una palabra binaria de 8 bits, con 0 correspondiente al negro y 255 a blanco. La resolución espacial de un píxel es de unos 80 m x 80m. Cada imagen contiene 2340 x 3380 de tales píxeles.

La base de datos es una subárea (pequeña) de una escena, que consta de 82 x 100 píxeles. Cada fila del set de datos se corresponde a una región de píxeles cuadrados 3x3 completamente contenida dentro de la subárea 82x100. Contiene los valoresen las cuatro bandas espectrales (convertidas en ASCII) para cada uno de los 9 píxeles en la región 3x3 y un número que indica la etiqueta de clasificación del píxel central.

La clase de cada píxel se codifica como un número, siendo estos:

1: suelo rojo
2: cultivo de algodón
3: suelo gris
4: suelo gris húmdeo
5: suelo con rastrojo de vegetación
6: mezcla de suelos
7: suelo gris muy húmedo
- N.B: No hay ejemplos con la clase 6 en este conjunto de datos (fueron removidos debido a dudas respecto a la validez de esta clase). Los datos se dan en orden aleatorio y se han eliminado ciertas líneas de datos para que no se pueda reconstruir la imagen original a partir de este conjunto de datos.

En cada línea de datos, los cuatro valores espectrales para el píxel superior izquierdo se dan primero, seguidos por los cuatro valores espectrales para el píxel medio superior y luego los del píxel superior derecho, y así sucesivamente de modo que los píxeles se leen en secuencia de izquierda a derecha y de arriba a abajo. Por lo tanto, los cuatro valores espectrales para los píxeles centrales están dados por los atributos 17,18,19 y 20. Si lo desea, puede usar solo estos cuatro atributos, ignorando los demás.

DATOS TOTALES → set de entrenamiento: 4435 set de prueba: 2000
CANTIDAD DE ATRIBUTOS → 36 (4 bandas espectrales x 9 píxeles en una región)
ATRIBUTOS → numéricos en el rango {0,255}
CLASES → 6 en total: {1,2,3,4,5,7}

Procedimiento

Empezamos tomando como válida la hipótesis del enunciado, la cual sugiere que el análisis puede conducirse tomando solo el píxel central. Extraemos los datos de entrenamiento y prueba para el píxel central a dos Data Frames: sat_trn_df y sat_tst_df respectivamente. Estos tendrán los valores para el píxel central en las cuatro bandas espectrales así como la clase correspondiente, resultando en Data Frames de 5 columnas en total.
Realizamos un breve análisis de visualización de estos datos. Este nos permite comprender su estructura y distribución, así como estimar la similitud entre los datos de prueba y entrenamiento.
Construimos los distintos modelos aplicables a un problema de clasificación multi-clase. Puntualmente: Árbol de Decisión, Random Forest, Vecinos Más Cercanos, Regresión Logística Multinomial, Análisis Discriminante Lineal y Cuadrático. Primero entrenamos estos modelos contra el set de entrenamiento y los evaluamos contra el mismo set de entrenamiento realizando Validación Cruzada con \(K=10\). En la medida en que exista una alta similitud entre los datos de entrenamiento y prueba, es justo asumir que el error de clasificación estimado de esta forma será semejante al obtenido al evaluar los modelos contra el set de prueba.
Para mantener el informe de una longitud razonable, y en base a análisis exploratorios sobre el dataset y lo visto en clase, se decidió no incluir en el mismo los modelos obtenidos mediante Boosting y Bagging, presentando únicamente Random Forest y Árboles de Decisión.
Evaluamos los modelos contra el set de prueba, derivando distintas métricas como el Error de Clasificación y Tablas de Confusión y realizamos un análisis comparativo de los modelos. También presentamos las Curvas ROC en formato One Vs. All.
Por último, repetimos todo el análisis precedente ahora tomando en cuenta todos los píxeles a modo de corroborar la hipótesis del enunciado y derivamos observaciones y conclusiones generales de este análisis.

Obtención de los Datos

Dado que el enunciado menciona que pueden utilizarse solo los datos espectrales correspondientes al píxel central, buscamos extraerlos, junto con su clasificación, a una sola Data Frame.

file_trn = "SAT_trn.txt"
data_trn <- read.table(file_trn, sep = "", header=F)
sat_trn_df <- data_trn      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)
file_tst = "SAT_tst.txt"
data_tst <- read.table(file_tst, sep = "", header=F)
sat_tst_df <- data_tst      %>% 
  dplyr::select(17:20,37)  %>% 
  rename(B1=V17,
         B2=V18,
         B3=V19,
         B4=V20,
         C=V37)

sat_trn_df$C[sat_trn_df$C == 7] <- 6
sat_tst_df$C[sat_tst_df$C == 7] <- 6
sat_trn_df_numeric <- data.frame(sat_trn_df)
sat_tst_df_numeric <- data.frame(sat_tst_df)
class_names <- c("rojo", "algodón", "gris","gris húmedo", "vegetación", "gris muy húmedo")
class_values <- c(1,2,3,4,5,6)
class_map <- setNames(class_names, class_values)
sat_trn_df$C <- class_map[sat_trn_df$C]
sat_tst_df$C <- class_map[sat_tst_df$C]

Visualización de los Datos

Contenidos del Set de Datos

B1	B2	B3	B4	C
92	112	118	85	gris
84	103	104	81	gris
84	99	104	78	gris
84	99	104	81	gris
76	99	104	81	gris
76	99	108	85	gris
80	112	118	88	gris
80	107	113	85	gris
76	91	104	74	gris húmedo
76	95	100	78	gris húmedo

Podemos comprobar que tenemos valores acotados entre 0 y 255 para cuatro bandas espectrales y una categoría de suelo asignada en cada fila.

Estructura y Distribución de los Datos

Observamos que los valores espectrales para los 3 tipos de suelo gris son similares en todas las bandas. Mientras tanto, el algodon es fácilmente distinguible, así como también los píxeles con vegetación o color rojo en menor medida, que parecieran agruparse por separado en el plano formado por las bandas B1 y B2 (variables V17 y V18)
En consecuencia, es justo esperar un error de clasificación mayor entre estas tres clases de grises. Conversamente, es justo esperar que el error de clasificación para el algodón y el suelo rojo sea considerablemente más bajo dados sus valores distintivos en las cuatro bandas.
Respecto a los espectros, las variables V17 y V18 son relativamente similares (R=0.8), relación que existe a su vez entre V19 y V20. Es probable entonces que un proceso de selección de variables lleve a mejores resultados a la hora de predecir, lo cual probaremos más adelante.

Comparación entre Set de Entrenamiento y Prueba

Antes de empezar, queríamos asegurar también que las proporciones de cada clase y de las propiedades se mantuvieran relativamente uniformes entre el test de entrenamiento y de testeo, ya que sino la comparación entre los resultados se vería comprometida, dificultando el análisis.

Cada linea representa una observación para las 4 bandas. La normalización de los valores no es necesaria ya que \(B_1:B_4\) estan en la misma escala. Podemos ver que la distribución de suelo en los datos de entrenamiento y prueba es prácticamente la misma.
En consecuencia, no va a ser necesario tener cuidados especiales a la hora de entrenar o comparar los modelos en base a la muestra de testeo. Es decir, la performance en la clasificación de una clase en particular no va a cambiar significativamente el resultado global de la precisión o error, ya que su proporción en la muestra se mantiene uniforme en ambos sets y los valores asociados a cada clase también son estables.

Construcción de Modelos

Árbol de Decisión

Empezamos con el uso de un árbol de decisión como primer modelo, utilizando la librería rpart en incluyendo las 4 bandas asociadas al pixel central. El árbol fue podado de forma tal de minimizar el error de clasificación estimado por Validación Cruzada con K=10.

Como comentario, para todos los modelos se utilizó un valor semilla fijo para obtener los mismos resultados en las distintas corridas, y asegurar que los datos de validación cruzada surgen de los mismos datos entre los modelos.

set.seed(20)
n <- nrow(sat_trn_df)
# Entrenamiento del Modelo con 10-fold CV  
sat_rpart <-  rpart(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method = "class", cp=-1,xval=10)
sat_rpart_summary <- printcp(sat_rpart)

## 
## Classification tree:
## rpart(formula = C ~ B1 + B2 + B3 + B4, data = sat_trn_df, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
## [1] B1 B2 B3 B4
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   2.6167e-01      0   1.00000 1.00000 0.0084779
## 2   2.5335e-01      1   0.73833 0.76896 0.0097636
## 3   1.2905e-01      2   0.48498 0.48498 0.0095487
## 4   6.3634e-02      3   0.35593 0.35742 0.0088020
## 5   1.7544e-02      4   0.29230 0.29379 0.0082400
## 6   1.4570e-02      5   0.27475 0.28189 0.0081181
## 7   6.6905e-03      7   0.24561 0.24234 0.0076694
## 8   5.9471e-03     10   0.22034 0.23253 0.0075467
## 9   5.6497e-03     11   0.21439 0.22629 0.0074660
## 10  4.7577e-03     13   0.20309 0.21885 0.0073673
## 11  4.6090e-03     14   0.19833 0.20934 0.0072364
## 12  3.2709e-03     16   0.18912 0.20636 0.0071945
## 13  1.6354e-03     17   0.18585 0.19774 0.0070698
## 14  1.3381e-03     19   0.18258 0.19625 0.0070479
## 15  1.1894e-03     21   0.17990 0.19536 0.0070346
## 16  1.0903e-03     22   0.17871 0.19477 0.0070258
## 17  8.9206e-04     25   0.17544 0.19715 0.0070610
## 18  6.9382e-04     31   0.17009 0.19715 0.0070610
## 19  5.9471e-04     34   0.16800 0.19715 0.0070610
## 20  4.9559e-04     48   0.15968 0.19685 0.0070567
## 21  4.4603e-04     51   0.15819 0.19387 0.0070124
## 22  3.9647e-04     55   0.15641 0.19328 0.0070035
## 23  2.9735e-04     58   0.15522 0.19328 0.0070035
## 24  1.4868e-04     61   0.15433 0.20012 0.0071047
## 25  8.4958e-05     63   0.15403 0.20101 0.0071176
## 26  7.4338e-05     70   0.15343 0.20071 0.0071133
## 27  0.0000e+00     74   0.15314 0.19952 0.0070960
## 28 -1.0000e+00    197   0.15314 0.19952 0.0070960

plotcp(sat_rpart)

t_nodes <- which.min(sat_rpart$cptable[,4])
cp<-sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"CP"]
CVErr_rpart = sat_rpart$frame[1,'dev'] / n * sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"xerror"]
poda_rpart<-prune(sat_rpart,cp=cp)
rpart.plot(poda_rpart, cex = 0.5, extra = 0)

Random Forest

La continuación del caso anterior dicho de alguna manera, el siguiente modelo fue Random Forest, un conjunto de árboles. Más adelante se verá el uso de este modelo con todas las covariables, donde su capacidad para reconocer las variables más importantes va a ser más significativo (al tener 36 variables en vez de 4)

set.seed(20)
sat_rf<- randomForest(C~., data= sat_trn_df, 
                     mtry=sqrt(ncol(sat_trn_df)-1), importance=TRUE)

importance(sat_rf)

##            1        2         3        4         5        6
## B1 227.87903 29.96135 362.16076 58.21353  50.21856 48.37936
## B2 140.86987 43.00905  53.51638 37.62347 112.03571 41.93234
## B3  30.04423 10.40657  24.77468 15.11574  28.56435 32.31841
## B4  57.76868 24.80617  34.26451 35.82601  52.16778 74.12139
##    MeanDecreaseAccuracy MeanDecreaseGini
## B1            186.14425        1186.7218
## B2            119.11363         914.1809
## B3             48.88156         425.2869
## B4             91.94379         853.1770

varImpPlot(sat_rf)

CVErr_rf <- unname(sat_rf$err.rate[nrow(sat_rf$err.rate),1])
CVErr_rf

## [1] 0.1492672

Más allá de eso, pareciera notarse una preferencia por las variables V17 V18 y V20 en los ranking de importancia.

Vecinos Más Cercanos

El siguiente modelo a evaluar es el de Vecinos Más Cercanos. En un primer momento se decidio analizar este modelo únicamente con los datos del pixel central debido a preocupaciones respecto a la maldición de la dimensionalidad. Sin embargo, más adelante veremos su performance con el dataset completo también.

El principal valor a ajustar es la cantidad de vecinos a considerar para la clasificación, el cual fue obtenido de forma tal de maximizar la precisión estimada a partir de Validación Cruzada (K=10).

set.seed(20)
sat_knn <- train(C ~ ., data = sat_trn_df, 
                   method = "knn",
                   preProcess = c("center","scale"),
                   trControl =  trainControl(method = "cv", number = 10),
                   tuneGrid = expand.grid(k = 1:100),
                   metric = "Accuracy")
results <- sat_knn$results
ggplot(results, aes(x = k, y = Accuracy)) + 
  geom_line(color = "steelblue", size = 1.2) + 
  xlab("Valor de K") + 
  ylab("Exactitud") + 
  ggtitle("Rendimiento para diferentes valores de K") +
  theme_bw() + 
  theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        axis.title.x = element_text(face = "bold", size = 12),
        axis.title.y = element_text(face = "bold", size = 12),
        axis.text = element_text(size = 10),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

best_k <- results$k[which.max(results$Accuracy)]
CVErr_KNN <- 1-unname(sat_knn$results[best_k,"Accuracy"])

Del gráfico de arriba se ve una rápida mejora en los resultados a partir de K=11, con el máximo obtenido para K=45.

Regresión Logística Multinomial

Siguiendo con otra familia de modelos completamente difrente, se procedió a utilizar un Regresión Logística Multinomial. Durante el curso sólo utilizamos la regresión logística para el caso binario, pero nos pareció interesante incluirla en este trabajo.

summary(sat_mlr)

## Call:
## nnet::multinom(formula = .outcome ~ ., data = dat, decay = param$decay)
## 
## Coefficients:
##   (Intercept)        B1          B2          B3         B4
## 2   -9.921154 0.5867093 -0.63317495  0.29943995 -0.1122932
## 3  -36.648978 0.7617752 -0.02673514  0.02795903 -0.2534971
## 4  -10.253548 0.5886442 -0.07999703 -0.02004337 -0.2854203
## 5    7.594031 0.4558313 -0.52173889  0.21832973 -0.2099831
## 6    4.599893 0.5678520 -0.18940106 -0.05394639 -0.2887664
## 
## Std. Errors:
##   (Intercept)         B1         B2         B3         B4
## 2   1.6977543 0.05040409 0.04477873 0.05284647 0.04508065
## 3   1.2832527 0.03391253 0.03743092 0.04631104 0.04434272
## 4   0.8182467 0.03437754 0.03467115 0.04409015 0.04040293
## 5   0.9108139 0.03161938 0.03489501 0.04077884 0.03601165
## 6   0.7841282 0.03491517 0.03359844 0.04276092 0.03778148
## 
## Residual Deviance: 3755.425 
## AIC: 3805.425

predicted_mlr <- predict(sat_mlr, newdata = sat_trn_df)
CVErr_mlr <- mean(predicted_mlr != sat_trn_df$C)

Análisis Discriminante Lineal

Luego, seguimos con un análisis discrimante lineal (LDA). Decidimos realizar un análisis de Best Subset para definir las mejores variables a incluir en el modelo (entre aquellas del pixel central).

set.seed(20)
sat_LDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="lda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)

##  `stepwise classification', using 10-fold cross-validated correctness rate of method lda'.

## 4435 observations of 4 variables in 6 classes; direction: both

## stop criterion: improvement less than 1%.

## correctness rate: 0.56054;  in: "B2";  variables (1): B2 
## correctness rate: 0.79776;  in: "B1";  variables (2): B2, B1 
## correctness rate: 0.82616;  in: "B3";  variables (3): B2, B1, B3 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        6.78

sat_LDA

## method      : lda 
## final model : C ~ B1 + B2 + B3
## <environment: 0x00000262a78aa0e0>
## 
## correctness rate = 0.8262

CVErr_LDA <- 1-unname(sat_LDA$result.pm["crossval.rate"])

De los resultados se ve que la mejor performance se obtuvo cuando se eliminó la variable B4 (V20), incluyendo únicamente B1-B3 (V17-V19). Las dos primeras variables coinciden con aquellas con mayor importancia acorde con el modelo de Random Forest, aunque la última se ve invertida.

Análisis Discriminante Cuadrático

Del mismo modo, se evaluó también el Análisis Discriminante Cuadrático, con el cual se espera obtener ligeramente mejores resultados que el caso anterior dada las características del dataset y la cantidad de datos disponibles.

Nuevamente se procedió con el método Best Subset para seleccionar las variables a incluir en el modelo.

set.seed(20)
sat_QDA <- stepclass(C ~ B1 + B2 + B3 + B4, data= sat_trn_df, method="qda", criterion = "CR", direction = "both", improvement <- 0.01, fold = 10)

##  `stepwise classification', using 10-fold cross-validated correctness rate of method qda'.

## 4435 observations of 4 variables in 6 classes; direction: both

## stop criterion: improvement less than 1%.

## correctness rate: 0.58107;  in: "B4";  variables (1): B4 
## correctness rate: 0.80496;  in: "B1";  variables (2): B4, B1 
## correctness rate: 0.84849;  in: "B2";  variables (3): B4, B1, B2 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00        7.36

sat_QDA

## method      : qda 
## final model : C ~ B1 + B2 + B4
## <environment: 0x0000026294058500>
## 
## correctness rate = 0.8485

CVErr_QDA <- 1-unname(sat_QDA$result.pm["crossval.rate"])

En este caso también se encontró que utilizar 3 variables era el óptimo, pero las variables elegidas (B1, B2, B4) coinciden con las reconocidas como más importantes según Random Forest.

Comparación del Error de 10-fold CV

Habiendo construido todos los modelos, a continuación se presentan todos los errores estimados por validación cruzada.

models <- c("AD", "RF", "KNN", "RLM","ADL","ADC")
CVErrs <- c(CVErr_rpart,CVErr_rf,CVErr_KNN,CVErr_mlr,CVErr_LDA,CVErr_QDA)

CVErr_df <- data.frame(Label = models, Value = CVErrs)

CVErr_plot <- ggplot(CVErr_df, aes(x = Label, y = Value)) + 
  geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Validación Cruzada (K=10)", title = "Comparación de Precisión de Modelos") +
  theme(plot.title = element_text(hjust = 0.5))

CVErr_plot

Se puede ver que no existe realmente mucha diferencia entre los modelos. Todos tienen una performance similar, con un 15% de error de clasifcación a nivel global aproximadamente. Los mejores resultados se obtuvieron utilizando Vecinos Más Cercanos, aunque la diferencia (absoluta) con sus competidores es menor al 1% en precisión.

Evaluación de los Modelos contra el Set de Testeo

A continuación, se muestran la tablas de confusión y la precisión de los modelos al utilizar el Set de Testeo. Se esperaría que los resultados sean comparables con aquellos obtenidos por validación cruzada, tal vez con algo menor de precisión ya que los datos nunca fueron utilizados en la construcción del modelo.

Se incluyen también las Gráficas ROC en formato One Vs. All, las cuales permiten evaluar qué tan separables son las clases en nuestros modelos. También podrían utilizarse para ajustar la metodología de clasificación si alguna misclasición se considerase más importante que otra. Para este trabajo se asumió que, dada la proporción de cada clase en la muestra y el tipo de problema (clasificación de suelo), no era importante optimizar los límites para elegir una clase por sobre otra. Nuestra única variable de interés era el error de clasificación a nivel global.

En nuestro caso, se puede ver que los Valores de AUC son muy altos (mayores a 0.99), lo que implica que los modelos hacen un buen trabajo separando la clases de suelo. En general se puede vaer que la curva celeste (gris húmedo) es la que menor área bajo la curva tiene, lo que era esperable por su similitud con los otros grises.

Ese mismo resultado puede verse en la matrices de confusión. Todas la matrices son análogas entre sí, con una muy buena performance para la mayoría de las clases pero con valores llamativamente bajos (no supera el 60% en el mejor de los casos) para el suelo gris húmedo, el cual es confundido la mayoría de las veces por suelo gris o suelo gris muy húmedo, lo cual es lógico.

Árbol de Decisión

sat_rpart_pred <- predict(sat_rpart,sat_tst_df,type="class")
utils.conf_matrix(sat_tst_df,sat_rpart_pred)

utils.OVAROC(sat_tst_df,sat_rpart)

CLErr_rpart <- mean(sat_rpart_pred != sat_tst_df$C)

Random Forest

Regresión Logística Multinomial

Vecinos Más Cercanos

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

Comparación del Error de Clasificación de los Modelos

En el gráfico abajo se puede ver la comparación de los errores de clasificación de todos los modelos. Los resultados fueron los esperados, manteniéndose las tendencias vistas en base a los errores estimados por validación cruzada, aunque ligeramente mayores.

Nuevamente, la performance entre los modelos fue similar entre ellos, con KNN proveyendo las mejores predicciones.

models <- c("AD", "RF", "KNN", "RLM","ADL","ADC")
CVErrs <- c(CLErr_rpart,CLErr_rf,CLErr_knn,CLErr_mlr,CLErr_LDA,CLErr_QDA)

CVErr_df <- data.frame(Label = models, Value = CVErrs)

CErr_plot <-  ggplot(CVErr_df, aes(x = Label, y = Value)) + 
  geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Clasificación", title = "Comparación de Precisión de Modelos") +
  theme(plot.title = element_text(hjust = 0.5))

CErr_plot

Análisis Para el Set de Datos Completo

Para comprobar (o rechazar) la idea de que la clasificación del pixel central puede obtenerse únicamente con los datos del propio pixel únicamente, se decidió repetir el proceso incluyendo la variables asociadas a los píxeles contiguos. Es probable que información del entorno permita a los modelos clasificar al pixel central mejor en casos más ambiguos.

A continuación, entonces, repetimos el análisis anterior, con los comentarios perinentes cuando haya alguna diferencia en el procedimiento o los resultados.

sat_trn_full_df <- data_trn %>% rename(C=V37)
sat_tst_full_df <- data_tst %>% rename(C=V37)
sat_trn_full_df$C[sat_trn_full_df$C == 7] <- 6
sat_tst_full_df$C[sat_tst_full_df$C == 7] <- 6
sat_trn_full_df$C <- as.factor(sat_trn_full_df$C)
sat_tst_full_df$C <- as.factor(sat_tst_full_df$C)

Árbol de Decisión

set.seed(20)
# Entrenamiento del Modelo con 10-fold CV  
sat_rpart <-  rpart(C ~ ., data= sat_trn_full_df, method = "class", cp=-1,xval=10)
sat_rpart_summary <- printcp(sat_rpart)

## 
## Classification tree:
## rpart(formula = C ~ ., data = sat_trn_full_df, method = "class", 
##     cp = -1, xval = 10)
## 
## Variables actually used in tree construction:
##  [1] V1  V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V2  V20 V21 V22 V23 V24 V25 V26
## [20] V27 V28 V29 V3  V30 V31 V32 V33 V34 V35 V36 V4  V5  V6  V7  V8  V9 
## 
## Root node error: 3363/4435 = 0.75829
## 
## n= 4435 
## 
##             CP nsplit rel error  xerror      xstd
## 1   0.26167113      0   1.00000 1.00000 0.0084779
## 2   0.25334523      1   0.73833 0.76896 0.0097636
## 3   0.12905144      2   0.48498 0.48498 0.0095487
## 4   0.06363366      3   0.35593 0.35742 0.0088020
## 5   0.01754386      4   0.29230 0.29497 0.0082518
## 6   0.01486768      5   0.27475 0.28338 0.0081337
## 7   0.01159679      6   0.25989 0.25453 0.0078153
## 8   0.00683913      7   0.24829 0.24591 0.0077128
## 9   0.00669045      9   0.23461 0.23877 0.0076253
## 10  0.00624442     13   0.20458 0.23699 0.0076030
## 11  0.00490633     15   0.19209 0.21588 0.0073269
## 12  0.00356824     17   0.18228 0.21053 0.0072531
## 13  0.00297354     18   0.17871 0.20607 0.0071903
## 14  0.00267618     19   0.17574 0.20517 0.0071776
## 15  0.00237883     21   0.17038 0.20369 0.0071563
## 16  0.00223015     22   0.16800 0.19923 0.0070916
## 17  0.00208147     24   0.16354 0.19923 0.0070916
## 18  0.00178412     26   0.15938 0.19744 0.0070654
## 19  0.00168500     27   0.15760 0.19685 0.0070567
## 20  0.00163544     30   0.15254 0.19685 0.0070567
## 21  0.00148677     32   0.14927 0.19506 0.0070302
## 22  0.00133809     36   0.14332 0.19387 0.0070124
## 23  0.00118941     40   0.13797 0.19447 0.0070213
## 24  0.00104074     50   0.12608 0.19387 0.0070124
## 25  0.00089206     55   0.12073 0.19120 0.0069721
## 26  0.00077312     57   0.11894 0.19060 0.0069631
## 27  0.00074338     62   0.11508 0.19090 0.0069676
## 28  0.00059471     64   0.11359 0.19090 0.0069676
## 29  0.00049559     66   0.11240 0.19179 0.0069811
## 30  0.00029735     69   0.11091 0.19269 0.0069946
## 31  0.00019824     74   0.10943 0.19804 0.0070742
## 32  0.00014868     77   0.10883 0.19893 0.0070873
## 33  0.00000000     83   0.10794 0.19982 0.0071003
## 34 -1.00000000    157   0.10794 0.19982 0.0071003

plotcp(sat_rpart)

which.min(sat_rpart$cptable[,4])

## 26 
## 26

cp<-sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"CP"]
CVErr_rpart = sat_rpart$frame[1,'dev'] / n * sat_rpart$cptable[which.min(sat_rpart$cptable[,"xerror"]),"xerror"]
poda_rpart<-prune(sat_rpart,cp=cp)
rpart.plot(poda_rpart, cex = 0.5, extra = 0)

En el resultado final que se peude ver arriba, se puede ver que el arbol de decisión obtenido mediante poda incluye numerosas variables fuera del pixel central, aunque aquellas del pixel central (V17-V20) tienden a “abrir el árbol”. Es decir, el modelo comienza la distinción con los datos del pixel central, y realiza las clasificaciones finales utilizando información “extra” de los píxeles contiguos.

Random Forest

##             1         2         3         4         5         6
## V1  14.347479  6.523180 12.785265  8.878950  9.215792 13.385304
## V2  25.080362  8.155885 10.329826 16.965498 13.801273 14.409123
## V3   7.750380  6.264034 10.674574 11.903256  8.868406  8.914603
## V4  11.671033  9.196722 16.190272 10.119210  9.599166 15.632320
## V5  20.327046  7.747617 14.732451  8.738162 12.260005 12.968500
## V6  22.024864  9.759857  8.967969 12.434673 15.591975 14.986900
## V7   8.193229  6.455881 11.596530 10.544200  7.867005  7.465667
## V8  16.440057  8.843152 10.938300  9.559178  9.546703 12.694114
## V9  27.245723  7.757990 11.797470  7.866552  9.557962 13.395249
## V10 19.668841 10.092615  7.702224  9.208108 14.255990 14.580319
## V11  8.183801  7.410958 12.179594 13.334420  5.674602 10.092250
## V12 13.869129  9.431361 12.276045 11.954095 10.082696 12.798514
## V13 20.251935  8.567868 16.168249  8.849309 11.890796 13.474792
## V14 22.215691 10.724722 10.223932 15.692282 15.305327 13.701134
## V15  7.862138  7.007398  9.906889 11.459271  7.941938 11.304583
## V16 16.505530 11.504687 13.949191  9.006578 11.792044 18.082491
## V17 25.253907 13.485706 19.712859 15.525689 20.078888 13.963230
## V18 21.396409 14.118719 14.385012 24.411052 24.468676 17.207413
## V19  8.560900  8.542848 10.437008 13.847183 10.239672 11.134022
## V20 19.027569 15.200911 11.888635 13.264401 11.666925 17.628338
## V21 29.431963  8.872517 17.675086 14.599193 15.461604 15.775755
## V22 20.051290 10.537004 11.771828 18.662235 19.376440 21.337874
## V23  9.749778  6.501513 12.096690 17.297252  8.046496  8.803671
## V24 15.582407  9.479060 14.049516 15.247126  9.821586 14.837700
## V25 23.697043  7.937527 16.986925  6.526053 10.340832 13.092715
## V26 23.867149  7.255891  9.029530 17.755496 10.805070 10.428707
## V27  7.894474  7.505076 12.865851  9.743226  8.409628  9.225300
## V28 14.727386 10.452824 13.915661  7.308338  9.939301 15.567830
## V29 22.313159  8.546988 15.018766 10.048277 10.401912 12.934963
## V30 17.037511  8.587731  8.721241 11.704029 13.527364 13.019518
## V31  8.545563  6.403730  8.403834 10.781731  9.028152  8.305877
## V32 16.176983  8.932304  8.691565 10.082833 11.848480 14.326251
## V33 24.911565  7.287095 15.507143  5.748705 10.478895 13.887847
## V34 18.162001  8.590430  9.937206  9.282485 16.150622 17.065096
## V35  8.769205  6.280038 12.924365 11.616752  8.099357 10.025106
## V36 16.279684  6.403133 13.129658  9.659609 11.145730 14.219610
##     MeanDecreaseAccuracy MeanDecreaseGini
## V1              19.84657         62.23829
## V2              32.28670         89.32196
## V3              18.44815         39.12481
## V4              20.91121         66.59387
## V5              22.37433        123.11917
## V6              31.22470        109.57090
## V7              15.57568         35.83360
## V8              21.46521         83.61104
## V9              27.50258        106.41434
## V10             29.38192         78.43511
## V11             18.88206         36.86583
## V12             23.55204         62.45378
## V13             21.89113        143.09708
## V14             30.72190        129.29556
## V15             18.39168         52.08460
## V16             25.45299        145.49467
## V17             27.32588        242.93899
## V18             33.44038        224.31144
## V19             16.77841         75.62461
## V20             23.97588        179.15180
## V21             30.97529        193.19072
## V22             33.09108        174.07645
## V23             21.39281         51.19399
## V24             23.28907        108.81224
## V25             26.31587        105.59507
## V26             30.82425         68.20054
## V27             20.38909         39.62732
## V28             21.72995         93.16480
## V29             24.21030        128.52837
## V30             27.27160         92.26471
## V31             17.24529         34.76598
## V32             20.34748         99.84895
## V33             29.07778        101.76329
## V34             29.70615        105.09069
## V35             21.32641         39.08788
## V36             24.29074         63.03664

La importancia de incluir todas las variables, y cuáles son de interés, es mas notorio en el modelo de Random Forest. Más adelante veremos el cambio en performance, pero más destacable en este momento es el gráfico de importancia de variables.

Viendo las gráficas se puede notar que, si bien el pixel central es determinante (V17, V18, V20, las mismas destacadas en los modelos anteriores y elegidas para QDA según Best Subset), otras variables juegan un rol importante en el modelo (V21, V22, V16).

Más adelante veremos que esto llevó a una reducción considerable en el error de clasificación.

Vecinos Más Cercanos

Siguiendo con KNN, en este caso se notó que el número de vecinos a incluir para maximizar la precisión se redujo considerablemente. Se puede ver un pico claro en 9 vecinos, mientras que antes se notaba una performance similar para un gran número de posibles valores de k.

Más adelante se verá también que KNN fue uno de los grandes beneficiarios de incluir más variables. Es probable que técnicas más eficientes de selección o reducción de variables puedan mejorar aún más su performance.

Regresión Logística Multinomial

## Call:
## nnet::multinom(formula = .outcome ~ ., data = dat, decay = param$decay)
## 
## Coefficients:
##   (Intercept)          V1          V2           V3           V4         V5
## 2  -10.041277 0.011633320 -0.08866011  0.048736082 -0.007836211 0.20621349
## 3  -20.378721 0.011140595 -0.06815353  0.027252249  0.039038401 0.12890553
## 4  -10.005767 0.037706249 -0.02832780  0.005182813  0.030101942 0.08771919
## 5    4.310293 0.015662639 -0.05979296  0.025596502  0.028523005 0.13158697
## 6    2.981347 0.004889374 -0.06399075 -0.020918896  0.065396066 0.15362516
##            V6         V7          V8          V9         V10         V11
## 2 -0.07046031 0.05849743 -0.09916763 -0.23035971  0.05487740 -0.01018132
## 3 -0.01023674 0.02334025 -0.10407334  0.01094643 -0.07388672  0.10356648
## 4 -0.04171813 0.07508490 -0.08062009 -0.03705449 -0.05554305  0.05291187
## 5 -0.03302893 0.02956949 -0.14203745 -0.04364464 -0.07723677  0.10014638
## 6  0.01751298 0.01713530 -0.13431599 -0.03026042 -0.04818043  0.04464404
##            V12        V13         V14          V15        V16          V17
## 2  0.047024992 0.11602469 -0.06530551  0.012422682 0.09757902  0.196116606
## 3 -0.053235834 0.12028057 -0.05625099 -0.013169008 0.06185362  0.021216840
## 4 -0.029957741 0.07085136 -0.05742333  0.002926769 0.03030945  0.042819710
## 5 -0.013539372 0.03250508 -0.04543212  0.028452811 0.05963729 -0.009482224
## 6  0.003462489 0.03033801 -0.04426977  0.005334445 0.04055678 -0.009438855
##            V18        V19        V20          V21        V22          V23
## 2 -0.170050812 0.08470339 -0.2117381 -0.025080470 0.09688088 -0.028190776
## 3 -0.024565151 0.09392422 -0.1519914  0.043435687 0.02603040 -0.035636777
## 4 -0.052699586 0.06156555 -0.1323629  0.094373655 0.03255762 -0.007513793
## 5 -0.088499090 0.10572731 -0.2332036  0.086437722 0.01248936 -0.032978130
## 6 -0.002384114 0.07188680 -0.1813734  0.009126433 0.02387681 -0.025139193
##          V24          V25         V26       V27         V28         V29
## 2 0.14146854 -0.104255928  0.01286808 0.1156963 -0.07636536 0.081602559
## 3 0.11258387  0.005678118 -0.06977875 0.1568363 -0.09402118 0.053380932
## 4 0.03780018 -0.085349494 -0.02746782 0.1701236 -0.13727885 0.035347957
## 5 0.08979900 -0.044821516 -0.03564611 0.1139013 -0.11675376 0.007046051
## 6 0.07970401  0.010761963 -0.04017290 0.0994432 -0.12406979 0.055388797
##            V30        V31          V32        V33         V34         V35
## 2 -0.111281819 0.04687024 -0.035408254 0.11878873 -0.07399343  0.01072422
## 3 -0.037891231 0.01333758 -0.043551988 0.08918813 -0.02290963 -0.01924823
## 4 -0.061578044 0.09632598 -0.073387776 0.12374137 -0.08545925  0.02117023
## 5 -0.006815922 0.03950480 -0.006685306 0.09955379 -0.09604781  0.02921400
## 6 -0.067981417 0.03953664 -0.012509769 0.07101754 -0.01432718  0.01335151
##             V36
## 2 -0.0410039083
## 3 -0.0007439179
## 4 -0.0184314429
## 5 -0.0439894601
## 6 -0.0664708960
## 
## Std. Errors:
##   (Intercept)         V1         V2         V3         V4         V5         V6
## 2  0.01236324 0.04366473 0.03165150 0.03499723 0.03335968 0.05167846 0.03677143
## 3  0.34938369 0.03234255 0.02519188 0.02874874 0.02940229 0.03930459 0.02912072
## 4  0.47599482 0.03283715 0.02525149 0.02873599 0.02899265 0.04009069 0.02926500
## 5  0.21208527 0.03591308 0.02666621 0.02943234 0.02900005 0.04281700 0.03056582
## 6  0.46297310 0.03106560 0.02392134 0.02651756 0.02719780 0.03737165 0.02767196
##           V7         V8         V9        V10        V11        V12        V13
## 2 0.03733509 0.03632063 0.04617490 0.03001791 0.03639816 0.03208830 0.04764122
## 3 0.02961524 0.03253631 0.03356620 0.02303724 0.02866407 0.02740302 0.03573847
## 4 0.02978146 0.03212779 0.03435290 0.02316587 0.02881193 0.02726241 0.03650879
## 5 0.03103002 0.03225148 0.03707527 0.02425580 0.02974641 0.02713716 0.03961988
## 6 0.02748290 0.03039313 0.03195555 0.02175804 0.02653014 0.02534765 0.03453444
##          V14        V15        V16        V17        V18        V19        V20
## 2 0.03403094 0.03590879 0.03675849 0.05359208 0.03706961 0.03792888 0.03849453
## 3 0.02676861 0.02912685 0.03248246 0.04044278 0.02918168 0.02968905 0.03469929
## 4 0.02668878 0.02931989 0.03221597 0.04089251 0.02916575 0.03003330 0.03408425
## 5 0.02835822 0.03012988 0.03222492 0.04407052 0.03128037 0.03112775 0.03408286
## 6 0.02570572 0.02700195 0.03036498 0.03850554 0.02797236 0.02750367 0.03269578
##          V21        V22        V23        V24        V25        V26        V27
## 2 0.04834975 0.03228657 0.03700916 0.03470027 0.04309672 0.02927003 0.03513068
## 3 0.03655183 0.02639977 0.02891477 0.03103616 0.03285555 0.02365707 0.02859243
## 4 0.03738730 0.02649046 0.02922960 0.03027532 0.03350740 0.02355047 0.02886155
## 5 0.04043290 0.02791366 0.03021617 0.03024462 0.03589307 0.02509726 0.02983479
## 6 0.03477958 0.02531165 0.02676384 0.02889859 0.03143473 0.02279398 0.02667284
##          V28        V29        V30        V31        V32        V33        V34
## 2 0.03391676 0.05116816 0.03648290 0.03758016 0.03769108 0.04256800 0.02926133
## 3 0.02951688 0.03896867 0.02922321 0.03015095 0.03424834 0.03124130 0.02349621
## 4 0.02929196 0.03953727 0.02916499 0.03031922 0.03362662 0.03277022 0.02353478
## 5 0.02966845 0.04199949 0.03090811 0.03129439 0.03337880 0.03614112 0.02492615
## 6 0.02750901 0.03670997 0.02787971 0.02790729 0.03212398 0.03059734 0.02249722
##          V35        V36
## 2 0.03600778 0.03289523
## 3 0.02876995 0.02939106
## 4 0.02880899 0.02859065
## 5 0.02978972 0.02839475
## 6 0.02666640 0.02697115
## 
## Residual Deviance: 4740.439 
## AIC: 5110.439

Análisis Discriminante Lineal

Tanto para QDA como LDA, se repitió el análisis mediante Best Subset. Inicialmente se trabajo con un Forward Stepwise para acelerar el proceso, pero luego se vio que el costo de utilizar Best Subset no era imposible en términos de tiempo.

Más allá de eso, los resultados obtenidos en este caso particular fueron idénticos entre ambos análisis. Para el caso de LDA, el mejor modelo se obtuvo incluyendo 3 variables: V17, V18 y V15. La última de estas variables esta fuera del pixel central, mientras que las dos primeras son aquellas que consistentemente fueron identificadas como las más importantes por todos los modelos.

## correctness rate: 0.56054;  in: "V18";  variables (1): V18 
## correctness rate: 0.79776;  in: "V17";  variables (2): V18, V17 
## correctness rate: 0.82797;  in: "V15";  variables (3): V18, V17, V15 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00       39.81

## method      : lda 
## final model : C ~ V15 + V17 + V18
## <environment: 0x0000026294f77218>
## 
## correctness rate = 0.828

Análisis Discriminante Cuadrático

El proceso se repitió para QDA. Lo interesante aquí es que se obtuvo exactamente el mismo modelo de antes! Con V17, V18 y V20. Los resultados son consistentes con el ranking de importancia de Random Forest, así como con los análisis realizadas en la sección anterior.

## correctness rate: 0.58107;  in: "V20";  variables (1): V20 
## correctness rate: 0.80496;  in: "V17";  variables (2): V20, V17 
## correctness rate: 0.84849;  in: "V18";  variables (3): V20, V17, V18 
## 
##  hr.elapsed min.elapsed sec.elapsed 
##        0.00        0.00       47.56

## method      : qda 
## final model : C ~ V17 + V18 + V20
## <environment: 0x0000026296cfcbe8>
## 
## correctness rate = 0.8485

Comparación del Error de 10-fold CV

Comparando la estimación de los errores obtenidos por validación cruzada, se destacan dos tendencias muy marcadas.

En general, los modelos de pocas variables (QDA, LDA, Regresión Logística) no vieron mejoras considerables en los resultados. Es más se mantienen en el orden de 15% de error. Esto era esperable, especialmente para QDA ya que “el mejor modelo” resultó ser el mismo independientemente del pool de variables con el que se empezó el Best Subset.

Lo que sí es notable es la enorme mejora de los resultados obtenidos mediante Vecinos Más Cercanos y Random Forest, los cuales redujeron a aproximadamente la mitad el error. Random Forest en particular, parece hacer un gran trabajo en reconocer internamente las variables más influyentes y fue un gran benefeciario de incluir todas las variables en su análisis.

Evaluación de los Modelos contra el Set de Testeo

Abajo se pueden ver las matrices de confusión para cada uno de los modelos. Si se mira detenidamente los resultados para KNN y Random Forest, lo más notable en comparación con los resultados anteriores es la gran mejora en la precisión para clasificar el suelo gris húmedo. Ese gran aumento de precisión en la clase más dificil (pasando, por ejemplo, de 39% a 57% en Random Forest) es una de las grandes razones de la mejora en precisión. Similares resultados se pueden ver en el caso de Vecinos Más Cercanos, y para el error en la clasificación de suelos con vegetación (con menores cambios, pero aún significativos).

Mientras tanto, el resto de los modelos mantuvo básicamente la misma performance en todas las clases, y no se vieron beneficiadas por la mayor información proporcionada por los píxeles contiguos.

Árbol de Decisión

Random Forest

Regresión Logística Multinomial

Vecinos Más Cercanos

Análisis Discriminante Lineal

Análisis Discriminante Cuadrático

Comparación del Error de Clasificación de los Modelos

models <- c("AD", "RF", "KNN", "RLM","ADL","ADC")
CVErrs <- c(CLErr_rpart,CLErr_rf,CLErr_knn,CLErr_mlr,CLErr_LDA,CLErr_QDA)

CVErr_df <- data.frame(Label = models, Value = CVErrs)

CErr_plot_full <- ggplot(CVErr_df, aes(x = Label, y = Value)) + 
  geom_bar(stat = "identity", color = "black", fill = "red") + labs(x = "Modelos", y = "Error de Clasificación", title = "Comparación de Precisión de Modelos") +
  theme(plot.title = element_text(hjust = 0.5))

CErr_plot_full

En el resumen de arriba se ve claramente como la performance de KNN y Random Forest es ampliamente superior para este dataset, con errores de aproximadamente 9% (contra 15% del resto de los modelos). A su vez, los resultados son similares a aquellos estimados por validación cruzada.

El mejor modelo resulto ser el de Random Forest.

Observaciones y Conclusiones Generales

Comparando el Error de Clasificación de los modelos considerando unicamente el pixel central y considerando todo el set de datos, concluímos que la hipótesis del enunciado no es válida, ya que se obtiene una reducción considerable en el error de clasificación al considerar todos los datos, así como una mayor variabilidad entre modelos.

CVErr_plot$labels$title <- "Error de CV (píxel central)"
CErr_plot$labels$title <-  "Error de Clasificación (píxel central)"
CVErr_plot_full$labels$title <- "Error de CV (completo)"
CErr_plot_full$labels$title <- "Error de Clasificación (completo)"

grid.arrange(CVErr_plot,
            CErr_plot,
            CVErr_plot_full,
            CErr_plot_full, 
            ncol = 2)

Considerando todos los datos, el mejor resultado se obtiene para Random Forest con un error de clasificación de aproximadamente el 8%, seguido de KNN con un error de clasificación de aproximadamente el 9%.
La mayor parte del error de clasificación viene dad por la dificultad de los modelos de distinguir la clase *gris húmedo de las otras categorías de grises, dato que se ve sustentado tanto análizando las matrices de confusión como las Curvas ROC One Vs. All.
Incluir los datos de píxeles aledaños ofereció una fuerte mejora de los resultados para KNN y Random Forest, el cual se explica principalmente por la mejora en la capacidad de distinguir la clase *gris húmedo, de 39% a 57% por ejemplo para Random Forest. En menor medida, la mejora en la clasificación de Vegetación también fue un factor importante.
Visto todo lo anterior, se selecciona como modelo el modelo de Random Forest incluyendo toda la información disponible para su construcción y evaluación.

TP Integrador - Aprendizaje Estadístico